Esse estudo foi realizado pelo grupo de estudos em ecologia espacial da UFBA, incluindo os seguintes pesquisadores:
Esse documento traz os resultados e códigos utilizados nas análises dos fatores que contribuem para a expansão dos casos de covid-19 no estado da Bahia.
Especificamente buscamos responder as seguintes perguntas nesse documento:
Como a taxa de crescimento tem variado ao longo do tempo na Bahia.
Quais fatores estão correlacionados à ocorrência de COVID-19 em um determinando município.
Quais fatores estão correlacionados ao número de casos de COVID-19 a cada 100 mil habitantes em um determinando município.
Quais fatores estão correlacionados ao número de dias até o primeiro e décimo caso de COVID-19 em um determinando município.
Apesar do número de casos terem subido, a taxa de crescimento do vírus na Bahia e em Salvador vem diminuindo, mas está aparentemente estabilizada no momento.
A presença de COVID-19 em um município parece ser influenciado pela presença de uma população grande e urbana. Municípios que estão próximos aos grandes aeroportos do estado em Salvador e Ilheus, também tem maior vulnerabilidade. Também é interessante notar que uma maior preciptação parece estar ligada a presença do vírus.
Em cidades que já tem a COVID-19, o número de casos por 100 mil habitantes parece aumentar principalmente com a proximidade com grandes aeroportos. O número de profissionais de saude também está positivamente correlacionado com o número de casos por 100k habitantes, indicando um potencial vies causado pela busca e disponibilidade de assitência medica. Os fatores ambientais também parecem ter um efeito sobre essa variável, indicando que o número de casos aumenta com maior precipitação e menor temperatura. Interessante é que a centralidade rodoviária parece ter um efeito negativo sobre o taxa de casos por 100 mil habitantes.
O número de dias até o primeiro e até o décimo caso é principalmente afetado pelo tamanho da população. A proporção de população rural parece diminuir o número de dias até a primeira contaminação também.
A versão do R utilizada foi:
R.version
## _
## platform x86_64-apple-darwin15.6.0
## arch x86_64
## os darwin15.6.0
## system x86_64, darwin15.6.0
## status
## major 3
## minor 6.2
## year 2019
## month 12
## day 12
## svn rev 77560
## language R
## version.string R version 3.6.2 (2019-12-12)
## nickname Dark and Stormy Night
Se for seguir o código para recriar as análises, antes de inciar, carregue e instale os seguintes pacotes.
library(coronabr) # pode baixar aqui: https://github.com/liibre/coronabr
library(tidyverse)
library(car)
library(randomForest)
library(rgdal) # load map
library(sp) # plot maps
library(plotly)
library(shiny)
Com o código abaixo podemos baixar os dados para todos os municipios Bahia. Para saber mais sobre as fontes do dados acesse o seguinte link: https://github.com/liibre/coronabr.
covid0 <- as_tibble(get_corona_br(uf = "BA"))
Pequenos ajustes na tabela:
covid <- covid0 %>%
filter(place_type == "city") %>%
mutate(city = factor(city, levels = unique(city)))
Dados por municipio:
mun_covid <- covid %>%
filter(date == date[1]) %>%
mutate(afetados = ifelse(confirmed > 0, 1, 0))
Estatísticas dos casos na Bahia:
mun_covid %>%
summarise(
"Casos totais" = sum(confirmed),
"Mortes totais" = sum(deaths),
"Número de municipios afetados" = sum(confirmed > 0)
)
Casos por município:
g <- mun_covid %>%
filter(confirmed > 0) %>%
mutate(city = as.character(city),
city = fct_reorder(city, -confirmed)) %>%
top_n(30, confirmed) %>%
ggplot(aes(x = city, y = confirmed)) +
geom_col() +
theme(axis.text.x = element_text(angle = 90)) +
scale_y_log10() +
xlab("") +
ylab("Casos confirmados de COVID-19") +
ggtitle("As 30 cidades mais afetadas na Bahia em número total")
ggplotly(g)
g <- mun_covid %>%
filter(confirmed > 0) %>%
mutate(city = as.character(city),
city = fct_reorder(city, -confirmed_per_100k_inhabitants)) %>%
top_n(30, confirmed_per_100k_inhabitants) %>%
ggplot(aes(x = city, y = confirmed_per_100k_inhabitants)) +
geom_col() +
theme(axis.text.x = element_text(angle = 90)) +
xlab("") +
ylab("Casos de COVID-19 por 100k hab") +
ggtitle("As 30 cidades mais afetadas na Bahia em casos por 100k hab")
ggplotly(g)
Número de casos no tempo.
covid_ba <- covid0 %>%
filter(place_type == "state")
g <- ggplot(covid_ba, aes(y = confirmed, x = date)) +
geom_point() +
geom_line() +
ylab("Casos confirmados") +
xlab("Data") +
ggtitle("Bahia")
ggplotly(g)
Aqui verificamos como a taxa de crescimento do vírus tem variado ao longo do tempo na Bahia.
rs <- covid_ba$confirmed[(nrow(covid_ba) - 1):1] / covid_ba$confirmed[nrow(covid_ba):2]
taxa <- tibble(data = rev(covid_ba$date[-nrow(covid_ba)]), rs)
g <- ggplot(taxa, aes(y = rs, x = data)) +
geom_smooth() +
geom_point() +
geom_line() +
scale_y_log10() +
ylab("Taxa de crescimento dos casos de COVID-19") +
xlab("Data") +
ggtitle("Bahia")
ggplotly(g)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Verifique o crescimento para cada município:
##
## Listening on http://127.0.0.1:4381
Aqui repetimos a análise para Salvador.
covid_sa <- covid0 %>%
filter(city == "Salvador")
rs <- covid_sa$confirmed[(nrow(covid_sa) - 1):1] / covid_sa$confirmed[nrow(covid_sa):2]
taxa <- tibble(data = rev(covid_sa$date[-nrow(covid_sa)]), rs)
g <- ggplot(taxa, aes(y = rs, x = data)) +
geom_smooth() +
geom_point() +
geom_line() +
scale_y_log10() +
ylab("Taxa de crescimento dos casos de COVID-19") +
xlab("Data") +
ggtitle("Salvador")
ggplotly(g)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'
Variação da taxa ao longo do tempo por municípios com mais de 20 casos confirmados. A variação é contada a partir do caso 10. Em azul, os muncipios onde a taxa vem diminuindo ao longo do tempo, em vermelho, os municipios onde a taxa vem aumentando.
cidades <- as.character(mun_covid$city[mun_covid$confirmed > 20])
n <- length(cidades)
rs_mun <- numeric(n)
for (i in 1:n) {
mun_i <- covid %>%
filter(city == cidades[i],
confirmed > 10)
if (nrow(mun_i) > 2) {
rs <-
mun_i$confirmed[(nrow(mun_i) - 1):1] / mun_i$confirmed[nrow(mun_i):2]
rs <- rs[!is.na(rs) & !is.infinite(rs)]
if (!all(rs == 1)) {
x <- 1:length(rs)
rs_mun[i] <- lm(rs ~ x)$coefficients[2]
}
}
}
taxa_mun <- tibble(cidades, taxa = rs_mun)
taxa_mun <- taxa_mun[rs_mun != 0, ]
g <- taxa_mun %>% mutate(cidades = fct_reorder(cidades, taxa)) %>%
na.omit() %>%
ggplot(aes(x = cidades, y = taxa, fill = taxa)) +
geom_col() +
theme(axis.text.x = element_text(angle = 90),
legend.position = "none") +
scale_fill_gradient2(low = "blue", mid = "gray", high = "red") +
xlab("") +
ylab("Variação da taxa")
ggplotly(g)
A primeira pergunta que buscamos identificar os fatores que afetam a ocorrência de COVID-19 em um determinando município. Buscamos fazer isso correlacionando diversos fatores sociais, econômicos e ambientais a duas variáveis respostas: uma binária, (1) separando municípios em infectados ou não infectados; e duas contínuas, (2) usando o número total de casos até o momento por 100 mil habitantes; e (3) identificando quanto tempo levou para o município ter X infecções por COVID-19.
No primeiro passo, carregamos os dados socio-economicos obtidos do IBGE, os dados de transporte obtidos do mapbiomas e dados do SUS.
ibge <- read_csv2(file("Data/new_ibge.csv", encoding = "UTF-8")) %>%
separate(cidade, c("Cidade", "Estado"), sep = "\\(") %>%
mutate(Estado = str_remove(Estado, "\\)")) %>%
filter(Estado == "BA") %>%
select(-X1, -X)
federal <- read_csv2("Data/federal_w_codes.csv") %>% select(-X1)
centralidade <- read_csv2("Data/new.dat.ba.csv") %>% select(-X1)
clima <- read_csv2("Data/climatic.br.csv") %>%
mutate(
precTotal = rowSums(.[3:7]),
tmean = apply(.[8:17], 1, mean)
) %>%
select(-X1)
meso <- read_csv("Data/meso.csv")
colnames(meso)[1] <- "mesoregiao"
aero <- read_csv2("Data/main.air.ba.csv") %>% select(-X1)
sus <- read_csv2("Data/new.data.sus.csv") %>% select(-X1)
dec <- read_csv("Data/decretos.csv") %>%
select(ibge, rod_fechada) %>%
filter(!duplicated(ibge))
Depois de carregados, juntamos as tabelas com todas essas informações.
munis <- left_join(ibge, centralidade, by = c("cod_ibge" = "ibgecode")) %>%
left_join(federal, by = c("cod_ibge" = "ibge")) %>%
left_join(mun_covid, by = c("cod_ibge" = "city_ibge_code")) %>%
left_join(clima, by = c("cod_ibge" = "ibge")) %>%
left_join(meso, by = c("cod_ibge" = "code")) %>%
left_join(aero, by = c("cod_ibge" = "ibge")) %>%
left_join(sus, by = c("cod_ibge" = "ibge")) %>%
left_join(dec, by = c("cod_ibge" = "ibge")) %>%
mutate(
dens.road = ifelse(is.na(dens.road), 0, dens.road),
afetados = ifelse(is.na(afetados), 0, afetados),
airport = ifelse(is.na(airport), "NO", airport),
confirmed = ifelse(is.na(confirmed), 0, confirmed),
rod_fechada = ifelse(is.na(rod_fechada), 0, rod_fechada),
leitos = ifelse(is.na(leitos), 0, leitos),
leitos = ifelse(is.na(leitos), 0, leitos),
profissionais.saude = ifelse(is.na(profissionais.saude), 0, profissionais.saude/total.pop),
confirmed_per_100k_inhabitants = ifelse(is.na(confirmed_per_100k_inhabitants), 0,
confirmed_per_100k_inhabitants
),
tam_pop_urb = total.pop * (1 - perc.rural)
)
Ao final obtemos a seguinte tabela:
munis
Com os dados de COVID-19 e os fatores economicos, sociais e ambientais, podemos explorar quais desses fatores afetam a COVID-19 no estado da Bahia. Iniciamos explorando quais desses fatores afetam a ocorrência de casos COVID-19 em um determinado município. Para tanto aplicamos um modelo logístico. Como muitos dos fatores estão correlacionados, escolhemos as variáveis com base nas hipóteses pensadas.
reg_log <- glm(afetados ~
airport + log(total.pop) + perc.rural +
log(eingen.cen.dist) + school.year +
perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
log(precTotal) + tmean,
data = munis,
family = binomial
)
Checar a inflação do modelo usando o valor de inflação de cada variável (VIF). A ideia é buscar manter o valor ao redor de 2 no máximo.
vif(reg_log)
## airport log(total.pop) perc.rural
## 1.250246 1.812368 1.669828
## log(eingen.cen.dist) school.year perc.with.wages
## 1.517274 2.215268 1.140648
## dist.min.ilh.ssa profissionais.saude log(precTotal)
## 1.352818 1.308366 1.206205
## tmean
## 1.197512
Para evitar possíveis problemas causados por desvios dos pressupostos do modelo logístico (homogeneidade e heterocedasticidade), recalculamos os valores de p por monte carlo, aleatorizando a variavel “afetados”.
resu <- summary(reg_log)
rept <- 1000
obs_z <- summary(reg_log)$coefficients[, 3]
obs <- coefficients(reg_log)
zs <- coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(zs) <- colnames(coefs) <- names(obs)
for (i in 1:rept) {
munis$rnd_afetados <- sample(munis$afetados)
reg_log <- glm(rnd_afetados ~
airport + log(total.pop)+ perc.rural +
log(eingen.cen.dist) + school.year +
perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
log(precTotal) + tmean,
data = munis,
family = binomial)
zs[i, ] <- summary(reg_log)$coefficients[, 3]
coefs[i, ] <- coefficients(reg_log)
}
for (j in 1:length(obs)) {
maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}
Aqui os resultados:
resu
##
## Call:
## glm(formula = afetados ~ airport + log(total.pop) + perc.rural +
## log(eingen.cen.dist) + school.year + perc.with.wages + dist.min.ilh.ssa +
## profissionais.saude + log(precTotal) + tmean, family = binomial,
## data = munis)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2973 -0.7414 -0.4628 0.7204 2.2153
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.082e+01 5.334e+00 -3.903 0.00200 **
## airportYES -9.861e-01 7.314e-01 -1.348 0.20180
## log(total.pop) 1.321e+00 2.548e-01 5.186 0.00200 **
## perc.rural -2.020e+00 8.429e-01 -2.396 0.01399 *
## log(eingen.cen.dist) 4.320e-03 1.023e-01 0.042 0.97103
## school.year -1.580e-01 2.888e-01 -0.547 0.57143
## perc.with.wages -1.131e+00 2.966e+00 -0.381 0.67532
## dist.min.ilh.ssa -3.573e-03 9.058e-04 -3.945 0.00200 **
## profissionais.saude 5.114e+01 5.247e+01 0.975 0.24975
## log(precTotal) 1.328e+00 4.668e-01 2.845 0.00599 **
## tmean 7.061e-02 1.173e-01 0.602 0.46753
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 508.37 on 404 degrees of freedom
## Residual deviance: 388.44 on 394 degrees of freedom
## (12 observations deleted due to missingness)
## AIC: 410.44
##
## Number of Fisher Scoring iterations: 5
g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>%
ggplot(aes(x = afetados ,
y = total.pop)) +
geom_violin() +
geom_jitter(alpha = .5) +
scale_y_log10() +
xlab("Município afetado") +
ylab("Tamanho população")
ggplotly(g)
g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>%
ggplot(aes(x = afetados ,
y = perc.rural)) +
geom_violin() +
geom_jitter(alpha = .5) +
xlab("Município afetado") +
ylab("Porcentagem população rural")
ggplotly(g)
g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>%
ggplot(aes(x = afetados,
y = precTotal)) +
geom_violin() +
geom_jitter(alpha = .5) +
scale_y_log10() +
xlab("Município afetado") +
ylab("Precipitação")
ggplotly(g)
g <- munis %>% mutate(afetados = ifelse(afetados == 1, "Sim", "Não")) %>%
ggplot(aes(x = afetados,
y = dist.min.ilh.ssa)) +
geom_violin() +
geom_jitter(alpha = .5) +
scale_y_log10() +
xlab("Município afetado") +
ylab("Menor distância para um grande aeroporto")
ggplotly(g)
A presença de COVID-19 em um município parece ser influenciado pela presença de uma população grande e urbana. Municípios que estão próximos aos grandes aeroportos do estado em Salvador e Ilheus, também tem maior vulnerabilidade. Também é interessante notar que uma maior preciptação parece estar ligada a presença do vírus.
No próximo passo, investigamos os fatores que afetam o número de casos confimardos para cada 100 mil habitantes.
reg_N <- lm(
log(confirmed_per_100k_inhabitants + 1) ~
airport + log(total.pop) + perc.rural +
log(eingen.cen.dist) +
perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
log(precTotal) + tmean,
data = filter(munis, confirmed > 0)
)
vif(reg_N)
## airport log(total.pop) perc.rural
## 1.426401 2.397254 2.031991
## log(eingen.cen.dist) perc.with.wages dist.min.ilh.ssa
## 1.735185 1.281817 1.315114
## profissionais.saude log(precTotal) tmean
## 1.277859 1.358267 1.269040
Recalculamos abaixo os valores de p usando uma simulação da Monte Carlo, onde o numero de casos foi aleatorizado e foi recalculada o número de casos por 100k habitantes.
resu <- summary(reg_N)
rept <- 1000
obs <- coef(reg_N)
coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(coefs) <- names(obs)
for (i in 1:rept) {
munis$rnd_cases <- (sample(munis$confirmed) / munis$total.pop) * 1000
while(sum(munis$airport[munis$rnd_cases > 0] == "YES") < 2) {
munis$rnd_cases <- (sample(munis$confirmed) / munis$total.pop) * 1000
}
reg_N <- lm(
log(rnd_cases + 1) ~
airport + log(total.pop) + perc.rural +
log(eingen.cen.dist) +
perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
log(precTotal) + tmean,
data = filter(munis, rnd_cases > 0)
)
coefs[i,] <- coef(reg_N)
}
for (j in 1:length(obs)) {
maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}
Aqui os resultados:
resu
##
## Call:
## lm(formula = log(confirmed_per_100k_inhabitants + 1) ~ airport +
## log(total.pop) + perc.rural + log(eingen.cen.dist) + perc.with.wages +
## dist.min.ilh.ssa + profissionais.saude + log(precTotal) +
## tmean, data = filter(munis, confirmed > 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4125 -0.5245 -0.0598 0.4748 1.8859
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.1101939 2.9462433 0.716 0.79321
## airportYES 0.3546829 0.3281526 1.081 0.15984
## log(total.pop) -0.1580030 0.1089460 -1.450 0.89111
## perc.rural -0.4113084 0.4714304 -0.872 0.20979
## log(eingen.cen.dist) -0.1067027 0.0615481 -1.734 0.00999 **
## perc.with.wages -0.0994676 1.8674938 -0.053 0.87512
## dist.min.ilh.ssa -0.0022055 0.0005062 -4.357 0.00200 **
## profissionais.saude 48.8380534 27.8826391 1.752 0.05395 .
## log(precTotal) 0.6773018 0.2601469 2.604 0.00400 **
## tmean -0.1109377 0.0671383 -1.652 0.02597 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7952 on 120 degrees of freedom
## Multiple R-squared: 0.3187, Adjusted R-squared: 0.2676
## F-statistic: 6.237 on 9 and 120 DF, p-value: 3.389e-07
g <- munis %>% filter(confirmed > 0) %>%
ggplot(aes(y = confirmed_per_100k_inhabitants,
x = eingen.cen.dist,
col = rod_fechada)) +
labs(colour = "Dias com a rodoviária fechada") +
scale_y_log10() +
scale_x_log10() +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Conectividade por transporte rodoviário")
ggplotly(g)
g <- munis %>% filter(confirmed > 0) %>%
ggplot(aes(y = confirmed_per_100k_inhabitants,
x = dist.min.ilh.ssa)) +
scale_y_log10() +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Menor distância para um grande aeroporto")
ggplotly(g)
g <- munis %>% filter(confirmed > 0) %>%
ggplot(aes(y = confirmed_per_100k_inhabitants,
x = precTotal)) +
scale_y_log10() +
scale_x_log10() +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Precipitação")
ggplotly(g)
g <- munis %>% filter(confirmed > 0) %>%
ggplot(aes(y = confirmed_per_100k_inhabitants,
x = tmean)) +
scale_y_log10() +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Temperatura")
ggplotly(g)
g <- munis %>% filter(confirmed > 0) %>%
ggplot(aes(y = confirmed_per_100k_inhabitants,
x = profissionais.saude)) +
scale_y_log10() +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Número de profissionais de saude")
ggplotly(g)
Em cidades que já tem a COVID-19, o número de casos por 100 mil habitantes parece aumentar principalmente com a proximidade com grandes aeroportos. O número de profissionais de saude também está positivamente correlacionado com o número de casos por 100k habitantes, indicando um potencial vies causado pela busca e disponibilidade de assitência medica. Os fatores ambientais também parecem ter um efeito sobre essa variável, indicando que o número de casos aumenta com maior precipitação e menor temperatura. Interessante é que a centralidade rodoviária parece ter um efeito negativo sobre o taxa de casos por 100 mil habitantes.
Na terceira linha de investigação, verificamos quais fatores afetam o númeri de dias até chegar ao caso, 1, 2, 3… N. Para tanto, usamos um modelo de linear com os mesmos fatores incluidos nos testes anteriores.
Antes de começar calculamos o número de dias até o primeiro caso por cidade, contados a paritir do dia 1 para o estado.
cidades <- na.omit(unique(covid$city_ibge_code))
n <- length(cidades)
baseline <- as.Date.character('2020-03-06')
dia1 <- numeric(n)
for (i in 1:n) {
covid_i <- covid %>%
filter(city_ibge_code == cidades[i], confirmed > 0)
dia1[i] <- as.numeric(min(covid_i$date) - baseline)
}
## Warning in min.default(structure(numeric(0), class = "Date"), na.rm = FALSE): no
## non-missing arguments to min; returning Inf
rem <- is.infinite(dia1)
dias <- tibble(cod_ibge = cidades, dia1)[!rem, ]
munis <- munis %>% left_join(dias)
## Joining, by = "cod_ibge"
## Warning: Column `cod_ibge` has different attributes on LHS and RHS of join
Modelo:
reg_dias <- lm(
dia1 ~
airport + log(total.pop) + perc.rural +
log(eingen.cen.dist) +
perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
log(precTotal) + tmean,
data = filter(munis, confirmed > 0
)
)
Verificando o efeito da colinearidade no modelo.
vif(reg_dias)
## airport log(total.pop) perc.rural
## 1.426401 2.397254 2.031991
## log(eingen.cen.dist) perc.with.wages dist.min.ilh.ssa
## 1.735185 1.281817 1.315114
## profissionais.saude log(precTotal) tmean
## 1.277859 1.358267 1.269040
Recalculamos abaixo os valores de p usando uma simulação da Monte Carlo, onde o numero de casos foi aleatorizado e foi recalculada o número de casos por 100k habitantes.
resu <- summary(reg_dias)
rept <- 1000
obs <- coef(reg_dias)
coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(coefs) <- names(obs)
for (i in 1:rept) {
munis_i <- filter(munis, confirmed > 0)
munis_i$rnd_cases <- sample(munis_i$dia1)
# while(sum(munis$airport[munis_i$rnd_cases > 0] == "YES") < 2) {
# munis_i$rnd_cases <- sample(munis_i$dia1)
# }
reg_dias <- lm(
rnd_cases ~
airport + log(total.pop) + perc.rural +
log(eingen.cen.dist) +
perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
log(precTotal) + tmean,
data = munis_i)
coefs[i,] <- coef(reg_dias)
}
for (j in 1:length(obs)) {
maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}
Aqui os resultados:
resu
##
## Call:
## lm(formula = dia1 ~ airport + log(total.pop) + perc.rural + log(eingen.cen.dist) +
## perc.with.wages + dist.min.ilh.ssa + profissionais.saude +
## log(precTotal) + tmean, data = filter(munis, confirmed >
## 0))
##
## Residuals:
## Min 1Q Median 3Q Max
## -22.5232 -6.2825 -0.9202 7.0412 21.4271
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.266e+02 3.761e+01 3.367 0.0360 *
## airportYES 3.043e+00 4.189e+00 0.726 0.5574
## log(total.pop) -5.866e+00 1.391e+00 -4.217 0.0020 **
## perc.rural 1.399e+01 6.019e+00 2.325 0.0599 .
## log(eingen.cen.dist) 1.347e+00 7.858e-01 1.715 0.1359
## perc.with.wages -2.051e+01 2.384e+01 -0.860 0.4056
## dist.min.ilh.ssa 4.076e-04 6.462e-03 0.063 0.9870
## profissionais.saude -3.890e+02 3.560e+02 -1.093 0.3996
## log(precTotal) -1.928e+00 3.321e+00 -0.580 0.5774
## tmean -1.035e-01 8.571e-01 -0.121 0.9211
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.15 on 120 degrees of freedom
## Multiple R-squared: 0.3488, Adjusted R-squared: 0.2999
## F-statistic: 7.141 on 9 and 120 DF, p-value: 3.046e-08
g <- munis %>% filter(confirmed > 0) %>%
ggplot(aes(y = dia1,
x = total.pop)) +
scale_x_log10() +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Tamanho da população")
ggplotly(g)
g <- munis %>% filter(confirmed > 0) %>%
ggplot(aes(y = dia1,
x = perc.rural)) +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Porcentagem da população em área rural")
ggplotly(g)
Antes de começar calculamos o número de dias até o décimo caso por cidade.
cidades <- na.omit(unique(covid$city_ibge_code))
n <- length(cidades)
baseline <- as.Date.character('2020-03-06')
dia10 <- numeric(n)
for (i in 1:n) {
covid_i <- covid %>%
filter(city_ibge_code == cidades[i], confirmed > 9)
dia10[i] <- as.numeric(min(covid_i$date) - baseline)
}
rem <- is.infinite(dia10)
dias <- tibble(cod_ibge = cidades, dia10)[!rem, ]
munis <- munis %>% left_join(dias)
## Joining, by = "cod_ibge"
Neste modelo precisaram ser retiradas algumas variáveis, devido a um aumento na inflação do modelo gerado pelo menor número de unidades amostrais.
munis_mod <- munis %>% filter(!is.na(dia10) & !is.infinite(dia10))
reg_dias <- lm(
dia10 ~
log(total.pop) + perc.rural +
dist.min.ilh.ssa +
log(precTotal) + tmean,
data = munis_mod
)
Verificando o efeito da colinearidade no modelo.
vif(reg_dias)
## log(total.pop) perc.rural dist.min.ilh.ssa log(precTotal)
## 2.077763 2.019820 1.836422 2.690450
## tmean
## 2.317169
Recalculamos abaixo os valores de p usando uma simulação da Monte Carlo, onde o numero de casos foi aleatorizado e foi recalculada o número de casos por 100k habitantes.
resu <- summary(reg_dias)
rept <- 1000
obs <- coef(reg_dias)
coefs <- matrix(ncol = length(obs), nrow = rept)
colnames(coefs) <- names(obs)
for (i in 1:rept) {
munis_i <- munis_mod
munis_i$rnd_cases <- sample(munis_i$dia10)
# while(sum(munis$airport[munis_i$rnd_cases > 0] == "YES") < 2) {
# munis_i$rnd_cases <- sample(munis_i$dia10)
# }
reg_dias <- lm(
rnd_cases ~
log(total.pop) + perc.rural +
dist.min.ilh.ssa +
log(precTotal) + tmean,
data = munis_i)
coefs[i,] <- coef(reg_dias)
}
for (j in 1:length(obs)) {
maior <- (sum(obs[j] >= coefs[, j]) + 1 ) / (rept + 1) * 2
menor <- (sum(obs[j] <= coefs[, j]) + 1) / (rept + 1) * 2
resu$coefficients[j, 4] <- ifelse(maior > menor, menor, maior)
}
Aqui os resultados:
resu
##
## Call:
## lm(formula = dia10 ~ log(total.pop) + perc.rural + dist.min.ilh.ssa +
## log(precTotal) + tmean, data = munis_mod)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.9123 -6.7459 0.8235 3.9126 15.5646
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 246.08333 108.52088 2.268 0.124
## log(total.pop) -7.95174 2.67545 -2.972 0.032 *
## perc.rural -3.19221 52.37697 -0.061 0.935
## dist.min.ilh.ssa -0.02800 0.04788 -0.585 0.713
## log(precTotal) -14.83370 18.57548 -0.799 0.615
## tmean -0.76395 3.11781 -0.245 0.875
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.554 on 10 degrees of freedom
## Multiple R-squared: 0.6138, Adjusted R-squared: 0.4207
## F-statistic: 3.179 on 5 and 10 DF, p-value: 0.05641
g <- munis_mod %>%
ggplot(aes(y = dia10,
x = total.pop)) +
scale_x_log10() +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Tamanho da população")
ggplotly(g)
g <- munis_mod %>%
ggplot(aes(y = dia10,
x = perc.rural)) +
geom_point() +
ylab("Número de casos por 100 mil habitantes") +
xlab("Porcentagem da população em área rural")
ggplotly(g)
Os resultados indicam que o tamanho da população é o principal fator determinando o número de dias até o primeiro e até o décimo caso. A proporção de população rural parece diminuir o número de dias até a primeira contaminação também.
Aqui analisamos